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We study the stability of spatially periodic, nonlinear Vlasov-Poisson equilibria as an eigenproblem 
in a Fourier-Hermite basis (in the space and velocity variables, respectively) of finite dimension, A'^. 
When the advection term in Vlasov equation is dominant, the convergence with A'^ of the eigenvalues 
is rather slow, limiting the applicability of the method. We use the method of spectral deformation 
introduced in [J. D. Crawford and P. D. Hislop, Ann. Phys. 189, 265 (1989)] to selectively damp 
the continuum of neutral modes associated with the advection term, thus accelerating convergence. 
We validate and benchmark the performance of our method by reproducing the kinetic dispersion 
relation results for linear (spatially homogeneous) equilibria. Finally, we study the stability of a pe- 
riodic Bernstein-Greene-Kruskal mode with multiple phase space vortices, compare our results with 
numerical simulations of the Vlasov-Poisson system and show that the initial unstable equilibrium 
may evolve to different asymptotic states depending on the way it was perturbed. 
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I. INTRODUCTION 

The stability of stationary nonlinear electrostatic 
waves, the so-called Bernstein-Greene-Kruskal (BGK) 
modes introduced in Ref. [1] , is a very old and basic prob- 
lem that is still of interest [2-4]. It appears, moreover, to 
be closely related to the saturation of stimulated Raman 
scattering (SRS), an issue that motivated the present 
work. Indeed, for physical conditions typical of those met 
in inertial confinement fusion, the electron plasma wave 
(EPW) that grows unstable due to SRS sees its ampli- 
tude change over space and time scales much larger than 
the Debye length and the plasma period, respectively [5] . 
In a sense, this EPW is therefore close to a BGK mode. 
Moreover, recent one-dimensional (1-D) Vlasov simula- 
tions of SRS [6] indicate that Raman reflectivity stopped 
increasing monotonically with time due to the growth of 
sidebands, resulting from a purely electrostatic instabil- 
ity similar to that introduced in Ref. [7]. More recently, 
the so-called vortex fusion instability [8], which we will 
further detail in this paper, was invoked in Ref. [9] to 
explain why Raman reflectivity stopped growing mono- 
tonically. 

In this paper, we strictly restrict to BGK equilibria, 
and describe a systematic and very efficient method to 
address their stability properties. Compared to a purely 
numerical approach consisting in integrating the Vlasov- 
Poisson system, our method not only very precisely pre- 
dicts the growth rate of the instability but also the func- 
tional form of the few fastest growing modes. This allows 
one to illuminate the physics behind the instability, to 
discern different pathways for subsequent evolution de- 
pending on the mode triggered and to devise viable con- 
trol strategies [10]. Although many approaches [2, 8, 11- 
13] have been developed to study the stability of non- 
linear electrostatic waves, to the best of our knowledge 
none has provided such a precise description of the un- 
stable modes by using a very general formalism and with 



very moderate computational cost. 

In order to determine the functional dependence of 
the unstable modes, an eigenproblem formulation is 
required. It is derived by linearizing the governing 
equations around any equilibrium distribution function 
fo{x,v). This leads to the following general formulation 

f^Af,. (1, 

where ^ is a linear operator which depends on /o, while 
fi{x,v,t) is an infinitesimal perturbation. The eigenval- 
ues of A determine the stability of the equilibrium char- 
acterized by /q. This equilibrium is unstable if some of 
the eigenvalues of A have a strictly positive real part, the 
largest of which being the growth rate of the instability. 

For spatially homogeneous equilibria, the eigenprob- 
lem (1) has been treated by Van Kampen [14] and 
Case [15]. For spatially inhomogeneous BGK equilibria, 
characterized by a non-vanishing electric field, we pro- 
pose here an approximate resolution of the eigenproblem 
(1) by making use of the Galerkin spectral method [16]. 
Hence, we expand the total distribution function over 
a finite set of global smooth functions which fulfill the 
boundary conditions of our problem, and which are more- 
over chosen to be orthogonal. Then, the operator A is 
approximated by a finite dimensional matrix. 

In many areas of physics, most notably fiuid mechan- 
ics, spectral methods have been particularly successful in 
solving eigenproblems of the form (1), with exponentially 
fast convergence of the result as a function of the num- 
ber of orthogonal functions retained in the expansion for 
a sufficiently smooth / (see Ref. [16], for example). This 
hallmark of spectral methods makes them attractive for 
the study of the stability of Vlasov-Poisson and Vlasov- 
Maxwell equilibria, especially in dimensions higher than 
one. However, a convincing application of spectral meth- 
ods to the problem of stability of BGK waves is still, to 
the best of our knowledge, lacking. Problems arise owing 
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to the very nature of the operator A. In particular the 
Unear, advection term in Vlasov equation [see Eq. (3a) 
below] is responsible for the transfer of energy to very fine 
velocity scales, a phenomenon known as velocity space 
filamentation (see, for example, Ref. [17]). When the 
contribution of the advection term is significant, eigen- 
functions of A are expected to involve high order velocity 
modes (fine velocity scales) leading to slow, power law 
convergence of the expansions on the orthogonal func- 
tions rather than the usual exponential convergence. An 
extreme example is the continuum of singular, neutral 
van Kampen modes whose approximation by means of 
smooth functions is clearly out of reach. Although we 
will not be interested in the continuum of neutral modes, 
it still has to be accounted for in the numerical approxi- 
mation of A and can interfere with the determination of 
the unstable modes of interest, particularly the weakly 
unstable ones. 

To ensure fast convergence of the eigenvalue calcula- 
tion, it is clear that the role of the neutral modes has 
to be undermined. To this end, we use the method of 
spectral deformation, originally developed for quantum 
mechanical problems and introduced in the study of the 
Vlasov-Poisson system by Crawford and Hislop [f8, 19]. 
The method introduces an operator 



B{e) = u{0)Au-^{e). 



(2) 



with U {9) non-unitary for complex 9. It can be 
proved [19] that the eigenvalues of A with nonzero real 
part remain unchanged under suitably chosen transfor- 
mations U{9), while the continuum of neutral modes is 
damped. Our central observation is that if the corre- 
sponding damping rate is chosen so that the continuum 
spectrum is well separated from the discrete eigenvalues 
of interest, then exponential convergence with the trun- 
cation order can be recovered. 

We choose to expand the distribution function in 
Fourier series (in space) and Hermite functions (in veloc- 
ity). This decomposition was first introduced by Grant 
and Feix [20] in the study of stability of spatially homoge- 
neous Vlasov-Poisson equilibria and recently generalized 
to the full Vlasov-Maxwell system and inhomogeneous 
equilibria by Camporeale et al. [21]. Even more recently 
Paskauskas and De Ninno [2] revisited the method from 
a nonlinear dynamics perspective, studying homogeneous 
equilibria of the Hamiltonian mean field model. We show 
through numerical examples that, with the introduction 
of spectral deformation, Fourier- Hermite expansions con- 
verge fast enough to be useful for the determination of 
unstable modes, even for strongly inhomogeneous equi- 
libria. 

The organization of this paper is as follows. We intro- 
duce the Vlasov-Poisson system and its linear approxima- 
tion around an inhomogeneous equilibrium in Sec. II. For 
simplicity, we restrict to one space and one velocity di- 
mension. Spectral deformation is introduced and applied 
to the Vlasov-Poisson system in Sec. III. In Sec. IV we 
derive the representation of B{9) in the Fourier-Hermite 



basis. We illustrate the treatment of Landau damped 
modes by our method in Sec. V A 1. In Sec. V A 2 we com- 
pare our results for a spatially homogeneous, bump-on- 
tail, distribution function against those obtained by nu- 
merically solving Landau's dispersion relation. We also 
use this test problem to illustrate some of the convergence 
issues that can arise as the wavelength of perturbations 
decreases and the advection term becomes more signifi- 
cant, and their resolution through spectral deformation. 
In Sec. V B we study a nonlinear example, namely a BGK 
mode with multiple phase-space vortices. Our results 
for the growth rate and the unstable modes agree with 
numerical simulations of the (nonlinear) Vlasov-Poisson 
system. At a qualitative level our calculations show that 
the collective modes trigger the vortex fusion observed 
in numerical simulations. We discuss our findings and 
the potential for optimization, as well as our future stud- 
ies based on this work, in Sec. VI. In the appendices we 
provide some technical details on the properties of the 
Hermite basis used here (Appendix A), the representa- 
tion of U{9) in Fourier-Hermite basis (Appendix B), the 
truncation of the representation of A (Appendix C) and 
the calculation of the Fourier-Hermite coefficients (Ap- 
pendix D). 



II. LINEARIZATION OF THE 
VLASOV-POISSON SYSTEM 

We consider the one-dimensional motion of electrons 
in an unmagnetized plasma, with immobile ions forming 
a neutralizing background. We restrict to situations that 
can be modeled using periodic boundary conditions in a 
spatial domain x £ [0, L]. The motion is described in 
terms of the Vlasov-Poisson system 



dt dx dv ' 



dE 
dx 



+ C30 



Jo 



fdv ~1 



dxE = 0, 



(3a) 



(3b) 



(3c) 



where velocity is normalized to the thermal one 



Vth 



a/2 



(eoTe/qlno 



1/2 



space to the Debye length Ad 



time to the inverse of the electron plasma 



,1/2 



and electric field to 



frequency ujpe = [qlno/eor 
Vfi^m^/ Xj^Qi,, where < is the electron charge. 

Let [fa{x, v), Eq{x)] be an equilibrium solution of (3a), 



^dfo 
dx 



E, 



dfo 
dv 



0. 



(4) 



The Vlasov equation is Galilean invariant and therefore 
any traveling wave solution can be reduced to an equilib- 
rium solution without loss of generality. 



3 



Substituting / = /o + /i , E = Eq + Ei in (3), where 
/i and El are infinitesimal perturbations, and accounting 
only for first order terms in /i, we get 



— ^-v^-(eo^+Ei^ 
dt dx \ dv dv 



dEi 
dx 



fi dv . 



dxEi^O, 



(5a) 



(5b) 



(5c) 



which describe the evolution of (/i,£^i) in the linear 
neighborhood of (/o, i?o)- 

Owing to periodic boundary conditions in the space 
variable, the distribution function and electric field may 
be expressed in terms of Fourier series. 



+ 00 

foix,v) = foiv)'^r{x)., 

r — — oo 



+ 00 



(6a) 



(6b) 



where ^r{x) = e"'''''^, /cq = 27r/L, and similar expansions 
hold for /i and Ei. Plugging system (6) into system 
(5) and using the standard Fourier basis orthogonality 
relations, we get 



d_ 
dt 



f^{v,t)^-ikkovf^{v,t) 



OO 



' 1 d 
fcn r dv 

1 d 

■ dv 



-T 

OO 



+ CXD 



dvfl{v,t) 



+ 00 



dvfSiv), (7a) 



mind that the electric field is determined self-consistently 
through Poisson's equation. 

Owing to the Hamiltonian structure of the Vlasov- 
Poisson system [22], eigenvalues = 7„ +iw„ of the real 
operator A come into quartets ±Ai, ±A* [see Fig. 1(a)]. 
Moreover, A characteristically has a continuum spectrum 
(Tc on the imaginary axis. 

If 7„ = Re(A„) > for some X„ G (T, a perturbation 
in the Ti'th eigenspace grows in modulus as e'''"', while 
its phase oscillates at frequency aj„ = Im(A„). Then, 
the norm of a generic perturbation having non-zero com- 
ponents along all eigendirections would grow asymptoti- 
cally in time as e™^'^'^'''"^*. In this case, we will say that 
the equilibrium is unstable. 

If no eigenvalue with strictly positive real part exists, 
Re(A„) = for all n, the eigenvalues form a continuoum 
which coincides with the imaginary axis and the corre- 
sponding eigenmodes are undamped (neutrally stable). 
These modes are singular (described by generalized func- 
tions or distributions) and do not represent physically 
observable modes of the system. However, their presence 
is connected to the collisionless damping of generic elec- 
tric field perturbations. Solving the initial value problem 
for small amplitude electrostatic waves in a Maxwellian 
plasma. Landau [23] famously showed that the electric 
field amplitude vanishes at an exponential rate. As shown 
by Van Kampen [14] and generalized to more general 
spatially homogeneous equilibria by Case [15], Landau 
damping can be understood as a destructive interference 
effect (known as phase mixing) of the neutral modes. 

In this paper, we will adopt the convention of sorting 
our eigenvalues by decreasing real part, so that for unsta- 
ble equilibria 71 corresponds to the largest growth rate 
and the eigenfunction ei is the fastest growing mode. 
Moreover, we will not use z = iA^ as is traditionally the 
case in plasma physics when writing solutions of (5). In- 
stead, we conform to the convention usually employed in 
the study of linear ordinary differential equations, which 
is more natural to a spectral discretization of the eigen- 
problem for A. As a result, the continuous part of the 
spectrum, CTc, lies on the imaginary axis, as illustrated in 
Fig. 1(a). 



if fc 0, 



(7b) 



III. SPECTRAL DEFORMATION 



where the prime in summations indicates that we omit 
the r = term, as we have incorporated (7b) into (7a) to 
eliminate the electric field. The restriction E^ = follows 
from condition (5c) on the electric field. Equation (7a) 
is of the form (1), 



dfi 
dt 



Ah 



where ,4 is a linear integro-differential operator that de- 
pends on /q. For the rest of this paper we will use the dis- 
tribution function alone to refer to solutions, keeping in 



As already stated in the introduction, in order to alle- 
viate the difficulty due to the continuum of neutral modes 
in the spectrum of A, we do not address the stability of 
a BGK equilibrium by solving the eigenproblem for A, 
but for the transformed operator 



with 



B{0) = u{e)Au~\e), 



(8) 



(9) 
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FIG. 1. (a) Spectrum of an unstable, spatially homogeneous 
Vlasov equilibrium: ±Ai and iA* are discrete eigenvalues, 
while (Jc is the continuous spectrum (gray, thickened for clar- 
ity), (b) Through spectral deformation, is shifted to the 
left half-plane, uncovering eigenvalues A2,A2 and A3, A3 cor- 
responding to damped modes [— Re(Ai) being the damping 
rate obtained through Landau's dielectric tensor formalism]. 
Since the Hamiltonian structure is destroyed by spectral de- 
formation, — Ai, — Ai disappear as CTc moves to their left (after 
Ref. [18]). 



Here, 



sgn(fc)6', k^O, 



0, 



k = 0. 



and h^{v) is the fcth Fourier mode of a given func- 
tion h{x,v). As shown by Hislop and Crawford in 
Refs. [18, 19], the main merit of this transformation is 
that for Im(^) < the continuous spectrum becomes 
damped, while the eigenvalues with Re(A) > remain 
unchanged [see Fig. 1(b) and Ref. [18] for a homogenous 
equilibrium]. When eigenvalues with Re(A) = of mul- 
tiplicity higher than one exist, we can talk of discrete 
eigenvalues "embedded" in the continuous spectrum (see 
Ref. [18]). These embedded discrete eigenvalues also re- 
main unchanged under spectral deformation. Hence, the 
stability issue of a BGK equilibrium may be equivalently 
addressed by calculating the eigenvalues of A or of B{9), 
except that these are much more easily and accurately es- 
timated for B{9), since the damped continuous spectrum 
is much easier to account for numerically. 

One may however wonder how Landau damping, which 
results from the phase mixing of the neutral modes, could 
be recovered using spectral deformation. As shown in 
Ref. [18], for a homogeneous equilibrium, when lm{9) < 
the neutral modes of A become damped by kg Im(6') . 
When this is larger in absolute value than the damp- 
ing rate 7„ obtained through Landau's dielectric tensor 



formalism, a new pair of discrete eigenvalues An, A* with 
I R-e(An)| = 7„ [and Im(A„) equal to the frequency pre- 
dicted by Landau] appears in the spectrum of B{9) [see 
Fig. 1(b)]. The complex conjugate eigenvalue corre- 
sponds to wavelength — /cq. Hence, Landau damping is in- 
deed recovered but, unlike for the original Vlasov-Poisson 
system, it now appears as being due to the damping of an 
eigenmode (corresponding to the largest negative eigen- 
value) for the dissipative dynamics represented by the 
operator B{9). This will be discussed in more detail in 
Sec. VAl. 

Let us now, as for the original Vlasov-Poisson system, 
write the eigenproblem for B{9) in Fourier space. Us- 
ing, (7a) and g{v) = U{9)f{v), one easily finds that the 
equation 



dgi 
dt 



= B{9)g^ 



(10) 



becomes, in Fourier space. 



d 

gi{v,t) = -ikko{v + ek)gi{v,t) 



dt 



+ 



ko ^ r dv 



-y' 



50 ''(" + 6k - Ok~r) 



°° 1 r> 

fco r ov 



dv gi{v,t) 
dvgUv). (11) 



In (11) the introduction of dissipation through spectral 
deformation with Im(^) < becomes apparent through 
the presence of 9k in the advection term. The invariance 
of the discrete spectrum with Re(A) > is not obvious 
from (11) but has been established in Ref. [19]. The 
unstable eigenmodes of the initial Vlasov-Poisson prob- 
lem can be recovered through the inverse transformation, 
f^{v,t) = U{-9u)g\{v,t), since dt{U-^g) = A{U-^g) 
and U (9) is time independent. 



IV. HERMITE EXPANSION 

In order to approximately solve the eigenproblem for 
B{9), we now write g'^{v) as a sum of orthonormal Her- 
mite functions. Such expansions have been first intro- 
duced in numerical studies of the Vlasov-Poisson system 

by Grant and Feix [20] and Armstrong [24] and present 
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various advantages. Hermite functions decay as 
at large v and therefore allow one to treat boundary 
conditions correctly without truncating the infinite in- 
terval. Moreover, convenient three term relations (see 
Appendix A) of the Hermite functions will result in 
an explicit, sparse representation of the operator B{9). 
Hermite functions are related to velocity moments, the 
few first of which are directly linked to physical observ- 
ables [20, 21]. As pointed out by Paskauskas and De 
Ninno [2], this allows one, at least in principle, to nat- 
urally separate thermal and filamentation scale effects. 
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As we will see in our numerical examples of Sec. V, fila- 
mentation scale often strikes back, rendering Hermite ex- 
pansions problematic, a shortcoming pointed out a long 
time ago [20] and overcome here by the introduction of 
spectral deformation. 

Here we consider the so-called asymmetrically weighted 
Hermite basis [25]. Denoting by ^'^(v) the basis func- 
tions and by ^'"(ti) the weight functions we have 



^„(t;) = C„e-" H^{v) , = CnHniv) , (12) 

where Hn{v) are Hermite polynomials and Cn = 
l/(7r^/^-\/2"ri!). More details on the properties of Her- 
mite polynomials used here can be found in Appendix A. 



We note the important orthonormality relation 

^+00 



(13) 



Our expansion of g'^{v,t) over the Hermite functions 
reads 



-t-00 

g''{v,t)=J2g'^m.{u), 

s=0 



(14) 



where v = v^u, with Vs an arbitrary velocity scale factor 
whose importance will be discussed at the end of this sec- 
tion. Plugging Eq. (14) into (11), multiplying by ^'^(u), 
integrating over u and using (A3), (A6) and the orthonor- 
mality of the Hermite basis we get, when j > 1, 



9i 



9i 



ITT-' ^/27_fc_,,,_l ,0 i^i/4 ggo ^1 



rl/4 



-5o 



r— — c)0 



fcn — ' r 

r— — oo n— 



(^'fc — 6'fe_,.) g\ , 



(15) 
(16) 



while 



dgl 



fcO 



dt 



-ikko 



9i +c^fc5i 



kO 



In Eq. (15) we have introduced the notations 



— > , if J < n, 



(17) 



[0, 

and g^^-'-iv) = g^^-^{v + 9k - 0k-r) = ft"{v + Ok) 
We can write (15) in tensorial form as 



dt 



where, 



kl 



V2 



5lm + Ok^On 



if j = „^ (18) while, when j > 1, 

if j > n, 



(19) and we have introduced 



(20) 



(21) 



-ikkoSki 



'j.m-^-l 



3 "I 



ka I 



<5omffo '"^ ^ ' if / ^ , 

if ; = , 



k-l.O 

9a 



9i) , if / 7^ fc and m < j — 1 
otherwise. 



(22) 
(23) 
(24) 



I 

Equations (20-24) provide the representation of the linear operator B{9) in the Fourier-Hermite basis. In 
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FIG. 2. Structure of 5*=^;™, for (a) 61 / and (b) 6 = 0, see 
(21) and (C2), respectively. We store B''-' im in a matrix as 
described in Appendix C. White space denotes zero elements. 



practice, we find it more convenient and efficient to 
compute the Fourier-Hermite coefficients /g^' of fo{x,v) 
rather than that of go{x, v). Then, (7q~''^ = /o " ^-^id we 

can compute gQ~^'''~^ = {U{Ok)fo)'^~^'''~^ as described in 
Appendix B. 

In computations, B''^ i,n has to be truncated to finite 
order by setting /^=^+^'^(t) = f''-^-+^{t) = for some 
cutoff values and iV^, see Appendix C. Physically, 
such a truncation holds provided that /o and fi can 
be well described by functions that do not oscillate too 
rapidly with space and velocity. The truncated matrix, 
shown in Fig. 2(a), is sparse and the computation of its 
first few eigenvalues with the largest real parts can be 
efficiently handled by iterative schemes, such as Arnoldi 
iteration [26]. For 6 = 0, A^^ im contains even fewer non- 
zero elements, see Appendix C and Fig. 2(b). 

The velocity factor Vs reflects the freedom to rescale 
the infinite domain in which Hermite functions are de- 
fined. This freedom has been exploited both in numeri- 
cal simulations [25] and eigenproblem calculations [2, 21]. 
Paskauskas and De Ninno use Vs to minimize the number 
of modes needed for a good approximation of the unper- 
turbed distribution function /q. Camporeale et ai, on 
the other hand, optimize Vg (and also allow for a shift 
of origin in u) iteratively so as to reduce the quadratic 
error in the eigenvectors /i, given an initial approxima- 
tion computed with non-optimal Vg ■ We have found that 
for the problems of interest to us, i.e. strongly inhomo- 
geneous equilibria, varying can be used to ensure fast 
convergence of the expansion of /q. However, in cases of 
power-law convergence of the eigenvalue computation, we 
have found that tuning Vg is of limited help as it would 
not yield exponentialy fast convergence. Therefore, in 
our numerical examples we will fix Vs to a value that pro- 
vides fast convergence of the expansion of /q and resort to 
spectral deformation to ensure exponential convergence 
rate of the eigenvalue calculation. 



V. NUMERICAL RESULTS 

A. Comparison with dielelectric tensor formalism 

For spatially homogeneous unperturbed distribution 
functions, the equilibrium electric field vanishes. The 
problem decouples in Fourier space and the growth (or 
damping) rate and frequency vary continuously with k 
and can be computed by Landau's dielelectric tensor for- 
malism [23] . We will test the validity of our expansions 
by comparing our results against the predictions of Lan- 
dau's theory for two examples. In the first example, /o is 
a Maxwellian, which lets us discuss the qualitative differ- 
ence in addressing Landau damping when diagonalizing 
A versus B{9). In the second example, we extensively 
compare the rate of convergence and the accuracy ob- 
tained with and without spectral deformation for an un- 
stable "bump-on-tail" distribution. 



1. Landau Damping 

Since we approximate the linear operator ^ by a finite 
dimensional matrix, the continuous, neutral spectrum is 
represented by a discrete set of eigenvalues, the separa- 
tion of which decreases as we increase Ny. As shown 
by Grant and Feix [20], Landau damping is then only 
approximately recovered when solving the initial value 
problem for the electric field. Indeed, if we consider a 
Maxwellian fo{v) and an initial small amplitude electric 
field perturbation of the form Ei(x,0) ecos/co^;, which 
is not an eigenmode of the problem, then the solution 
of the initial value problem in the linear approximation 
yields for the first Fourier component. El, of the electric 
field [20] 

E\(t) = El{Q)Y^{a-^),,a,,e^-^\ (25) 

where the w^-'s are the eigenvalues of A^^ [given by 
(C2)], aij the matrix of its column eigenvectors, and 
{a~^)ij the matrix of row vectors of the dual basis. 
Hence, Landau damping of the EPW occurs as the conse- 
quence of the destructive interference effect of the neutral 
eigenmodes, but only over a finite time, because the sum 
(25) is finite. This is illustrated in Figure 3(a) plotting 
the time evolution of the EPW amplitude as predicted 
by (25), for fcg = 0.5. The electrostatic wave is indeed 
initially damped and the damping rate and frequency of 
oscillation can be determined from the slope and local 
maxima distance in Fig. 3(a), 7 ~ —0.1534, w ~ 1.416. 
As shown in Fig. 3(b) and (c), agreement with the re- 
sults from Landau's analysis is excellent. However, after 
t ~ 24 the electric field amplitude grows again to finite 
magnitude, a fact known as recurrence. The recurrence 
time increases as we decrease the scale factor Vs or we 
increase Ny (see Ref. [25]) indicating that recurrence is 
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FIG. 3. (Color online) (a) Landau damping of an EPW as prescribed by Eq. (25), with = 51, Vs — \/2 (blue dotted line), 
spectral deformation with 9 — — 2.5i, A*',, = 51, Vs — \/2 (red solid line) and Landau's prediction for the damping of the envelope 
(black solid line). The recurrence effect is clearly seen in the blue dotted line for lm{9) — 0. (b) Landau damping rate and (c) 
frequency of oscillation of an EPW as predicted by Landau's analysis (black solid line) and Fourier-Hermite expansion with 
= 51, V, = ^/2, Im(6l) = (blue dots). 



ultimately related to the approximation of the continu- 
ous spectrum Uc by a finite set of eigenvalues (see also 
Ref. [27]). 

With spectral deformation, 9 — — 2.5i, the physics 
of the eigenproblem becomes dissipative: the Landau 
damped mode appears in the spectrum as a true eigen- 
mode with Ai ~ —0.1534 + 1.416 i, accompanied by 
and recurrence is absent. At the same time some of the 
collisionless physics is sacrificed, namely phenomena that 
depend on the reversible structure of Vlasov equation 
and the creation of fine velocity scales in the distribution 
function, such as plasma echos (See Ref. [28], for exam- 
ple). We do not plot the results obtained with spec- 
tral deformation in Fig. 3(b) and (c) since, within the 
resolution of this plot, they perfectly overlap with the 
Im(6l) =^ rcsuhs. 

2. Bump-on-tail instability 

Let us now compare the results of our method against 
those obtained from the numerical resolution of Landau's 
dielectric tensor formalism for a bump-on-tail distribu- 
tion of the form 

/(«) = ^e--^^ + ^e-(— (26) 

with Up = 0.9, rib = 0.1 and Vb = 5, see Fig. 4(a). 

The probability distribution function (26) can be effec- 
tively approximated with Hermite functions, as indicated 
by the exponential decay of the magnitude of Hermite 
coefficients f^-' in Fig. 4(a), which fall bellow roundoff 
foi' j ^ 140 (with Vs = 2.2). Note that both axes are 
logarithmic in the inset of Fig. 4(a) and exponential con- 
vergence corresponds to an envelope of the /q"* with an 
ever-increasing negative slope. For Ny = 60, the residual 
in /q"' is of the order of 10"''. In Fig. 4(b) we compare our 
results for 7 to numerical values obtained through Lan- 
dau's analysis. Without spectral deformation, lm(0) = 
and Ny = 59 and 60, agreement is good for small wave- 
lengths, but deteriorates as fcp is increased. In particu- 
lar, the mere addition of one Hermite term in the series 



changes the results significantly, and for Ny odd [20] the 
unstable feg range is narrower than predicted by Landau's 
analysis. Beyond the cutoff wavelength, the damping rate 
would have to be determined as in Sec. V A 1, but we will 
not get into this, as we are primarily interested in grow- 
ing modes. 

To study how the accuracy of our calculations is af- 
fected by Ny we fix /cq — 0-36 and vary the number 
of Hermite polynomials up to Ny = 200, see Fig. 4(c). 
As a measure of convergence we plot in Fig. 5(a) \jn^ — 
lN^-i\/lN^ , i.e. the relative change in 7 with the addition 
of a new term in the Hermite series, while in Fig. 5(b) we 
compare 7jv„ with the exact value of 7 derived from Lan- 
dau's analysis. Without spectral deformation we observe 
slow, power-law convergence. Beyond Ny — 140 the ex- 
pansion of /o only adds numerical noise to the eigenvalue 
computation; however the relative change in eigenvalues 
is of the order of 10~^, with pronounced oscillations with 
odd and even order truncation in Ny, see also Fig. 4(c). 

This even-odd order oscillation behavior has been ob- 
served by other authors in similar [20] or more general 
settings [21]. While some (see Ref. [2]) overcome such 
difficulties by averaging the eigenvalues computed over 
different values of Ny , such an approach does not warrant 
convergence and cannot justify the choice of an expansion 
over Hermite polynomials, as opposed to an estimation 
with a low order method, for instance with finite differ- 
ences. 

Equation (7) shows that the relative importance of the 
advection term, responsible for the poor convergence of 
the method, increases with feg. For relatively large feg's, 
spectral deformation is called for, and with 9 = — 2i and 
Ny = 60, agreement with Landau's analysis is recovered 
for all fco's in Fig. 4(b). The unstable range in fcg is 
now accuratelly retrieved, and the stability threshold is 
crossed smoothly, since damped modes are represented 
as true eigemodes of B{9) rather than through the inter- 
ference of neutral modes. 

A closer look at the convergence rate for different 
values of lm{9) in Fig. 5(a) reveals that for small val- 
ues of — lm(0) the convergence still obeys a power-law. 
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FIG. 4. (Color online) (a) Bump-on-tail distribution (26) with rip = 0.9, rib = 0.1, vt = 5 and the magnitude of its Hermite 
expansion coefBcients |/q-' | in the inset, (b) Growth rate 7 as a function of ko as predicted by Landau's method (red solid 
line) and by Hermite expansion, with and without spectral deformation (black lines), (c) Variation of the growth rate 7 with 
increasing A'^„, compared to the exact value for fco — 0.36 and — 0, — 0.5i, — 2i, — 3i. 




FIG. 5. (Color online) (a), relative change of 7]v„ with A'^„ and (b), relative difference to the exact value from Landau's 
analysis, for different values of lm{9). (c) Magnitude of Hermite expansion coefficients \fl-'\ of the eigenmode as computed with 
lm{9) = 0, -2. 



yet steeper than the one for lm{6) — 0. For large 
enough values of — lm(0), convergence becomes exponen- 
tial and there appears to be no practical advantage in 
further increase of — Ini(f?), since the convergence rate 
remains practically the same. From the results plotted 
in Fig. 4(c), we could say that the eigenvalue compu- 
tation follows an overdamped oscillation pattern, with 
very fast relaxation towards the exact value. As the ex- 
pansion coefScicnts /q"' fall below roundoff at j ~ 140, 
further precision gain with an increase in Ny ceases, see 
Fig. 5(a) and (b). 

It is interesting at this point to examine how rapidly 
the expansion coefficients fl'' of the computed eigenfunc- 
tion fall off. With either Im(6') = or Im(6i) < 0, coef- 
ficients fl'' fall off as a power-low, rather than exponen- 
tially as for /q-*, see Fig. 5(c). This should have been 
anticipated: owing to the action of the advection term 
the eigenfunctions span all velocity scales, up to the fil- 
amention scale. We therefore need non-vanishing contri- 
butions from higher order Hermite functions to capture 
such thin scale effects in the eigenfunction. The pay- 
back of using spectral deformation is a (relatively) ac- 
curate computation of high j spectral coefficients. On 



the contrary, with Im(6') = higher j spectral coeffi- 
cients are subject to strong even-odd order oscillations, 
see Fig. 5(c). 

B. Stability of BGK waves 

Bernstein-Greene-Kruskal (BGK) modes [1] are sta- 
tionary, nonlinear electrostatic waves accounting for the 
presence of both trapped and untrapped particles. They 
have attracted a lot of interest because of their resem- 
blance with the saturated state of nonlinear processes in 
plasmas. We will study the stability of a particular BGK 
equilibrium introduced by Ghizzo et al. [8] to model pro- 
cesses that involve vortex fusion. Although no claim can 
be made regarding the generality of the shape of /o pre- 
scribed in Ref. [8] , it serves well to demonstrate how our 
stability calculation can be used to predict vortex fusion. 

The chosen equilibrium distribution function is 

where H{x,v) — + $(2;) is the total energy, /i < 1 
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X X 

(c) (d) 

FIG. 6. (Color online) (a) Two spatial periods of the distribution function (27) with n = 0.92, ^ = 0.90 (solid lines indicate 
isocontours) , (b) the unstable eigenfunction of the 2-cell equilibrium, determined using A^^, = 48, Nv = 110, Vs — 1.8, = — i, 

(c) linear evolution in the neighborhood of the BGK mode, showing the holes approaching each other and triggering fusion, 

(d) late time {t ~ 1600) state of the initially perturbed equilibrium obtained with the Eulerian Vlasov code VADOR [29]. 



a parameter that controls inhomogeneity {fi — 1 corre- 
sponds to the homogeneous case), and ^ < 1 a parame- 
ter that controls the depth of the distribution function's 
"depression" or "hole" at x — (see Ref. [8] for details). 
Using Poisson equation, one then finds that the potential 
<i>(x) solves 

*"M = -.^^|^e-« + l. (28) 

In order to allow for subharmonic perturbations, 
Ghizzo et al. study in Ref. [8] the stability of the equi- 
librium (27) for a physical system whose periodicity is 
TV times the period, A, of the BGK mode. This sys- 
tem may be viewed as an "A^-cell replica" of the basic 
cell. Using the marginal stability analysis developped 
in [12], Ghizzo et al. show that when N > 2 the equi- 
librium is unstable. This theoretical prediction is then 
tested against long-time numerical integrations of the 
Vlasov-Poisson system, with an initial state consisting 
of the equilibrium (27) perturbed by a small amplitude 
monochromatic wave whose wavelength is L = NA. For 
TV = 1, the numerical results suggest that the equilibrium 
is stable. For N > 2, instability is clearly demonstrated 
in Ref. [8], as the perturbed equilibrium evolves towards 
a different final state through hole merging. 



Here, we compute the unstable modes and their 
growth rates for the BGK equilibrium (27) using our 
Galerkin projection method. Unless otherwise noted, for 
all Fourier-Hermite calculations that follow, we use = 
24 points/cell, Vg = 1.8, Ny in the range 40 - 110, 
and = 01 — i. Since growth rates are not com- 
puted in Ref. [8] , we compare our results with the growth 
rates obtained numerically from the resolution of the 
Vlasov-Possion system with the ID Eulerian Vlasov code 
VADOR [29]. 

Following Ref. [8], we solve Eq. (28) numerically for 
fi = 0.92, ^ = 0.90. The spatial period of the solution is 
fixed to A = 14.7106 by specifying the initial condition 
$(0) = $'(0) = 0. Phase space portraits of the distribu- 
tion function are characterized by one vortex or "hole" 
for each spatial period, as plotted in Fig. 6. Fourier- 
Hermite expansion coefficients of /o fall off exponentially, 
see Fig. 7(a). 

For N = 2, the Fourier-Hermite method converges 
rapidly to the growth rate 7 = 0.05562 ± 0.00001 when 
we use spectral deformation with 6 = — i, with the rel- 
ative error falling bellow 10~^ for iV^, > 90, see Fig. 8. 
The 6 — computation on the other hand suffers from 
large amplitude oscillations between even and odd order 
in Ny and the convergence rate is slow. 
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FIG. 7. (Color online) (a) Magnitude of expansion coeffi- 
cients I/q-'I of the BGK mode distribution function (27) for 
N = 2, computed with Vs = 1.8, N^: = 48. (b) Magnitude 
of expansion coefficients l/f"*! of the corresponding unstable 
eigenfunction. Only the leading Fourier modes are shown. 

The unstable eigenfunction is plotted in Fig. 6(b). In 
Fig. 6(c) we present the predicted linear evolution in the 
neighborhood of the BGK mode, showing that the in- 
stability causes the holes to approach each other. This 
triggers the eventual merging of the two holes, illustrated 
in Fig. 6(d), under the nonlinear dynamics. 

The distribution function corresponding to the unsta- 
ble mode determined by using the Fourier-Hermite ex- 
pansion with 9 = and for A^^, = 109 and 110 is shown 
in Fig. 9(a). We can see that the even/odd oscillations 
in Fig. 8(a) translate into rather large deviations in the 
shape of the eigenf unctions. 

Our results for 7 and the unstable mode, denoted 
as ei, are compared against those obtained through 
Vlasov simulations. We initialize the code VADOR 
with a perturbed distribution function of the form 
fo{x,v) [1 + acos(A;oa; -I- 7r/4)], with a = 10~^, and fol- 
low the evolution of the system in the neighborhood of /o 
until fi{x,v,t) = f{x,v,t) — fo{x,v) assumes a constant 
shape, plotted in Fig. 9(b), and only changes in norm. We 
estimate a growth rate 7„„m = 0.0556 ± 0.0002 from the 



rate of change of the distance \\f{x,v,t) — fo{x,v)\\2- 
Note that the resolution used in VADOR, = 1920, 
rijj = 1000, is much higher than what we needed in the 
Galerkin method, while the precision is lower, due to the 
error in graphically estimating the growth rate of the per- 
turbations. We compare the distribution functions corre- 
sponding to the predicted unstable mode with 9 = — i, to 
the one determined by numerical integration of a small 
sinusoidal perturbation with code VADOR in Fig. 9(b). 
We observe that the two profiles agree rather well, ex- 
cept for some small-scale structure not captured by our 
method. However, the resolution is too different in the 
two methods for a direct comparison at this scale to be 
meaningfull. 

In Fig. 7(b), we observe that the high-order Hermite 
function components of the unstable eigenmode do not 
fall off rapidly, indicating that the eigenf unctions we ap- 
proximate involve fine velocity scales. The utility of the 
combination of spectral deformation and Galerkin projec- 
tion is that it allows the accurate representation of the 
thermal scale effects of the instability, independently of 
filamentation scale effects. By contrast, not using spec- 
tral deformation {9 = 0) results in a much coarser ap- 
proximation of the eigenmodes, with large fluctuations 
between odd and even N^, see Fig. 9(a). However, both 
9 = —i and 9 = calculations provide very accurate 
approximation of the electric field, Fig. 9(c), as the dif- 
ferences in the distribution function are smoothed out 
when one considers its lower moments. 

The 3-cell system has a slightly smaller positive eigen- 
value 7 = 0.04856 ± 0.00001, with an eigenmode that 
leads to two-hole fusion, as shown in Fig. 10. The nu- 
merical simulations in Ref. [8] show that this is indeed 
the case, with the third hole subsequently merged with 
the other two, leading to an asymptotic one-hole state. 
We will not present here any detailed convergence study 
or comparisons with numerical simulations, as the results 
are qualitatively similar to the 2-cell case. 

For the 4-cell state a new possibility arises, since we 
can think of it as two 2-cell systems stacked together. 
Thus, the 2-cell, k = ko/2, unstable mode, denoted as 
ei, is still present, with eigenvalue 71 = 0.05562. There 
is yet another unstable mode, 82, with smaller eigen- 
value 72 = 0.04009, which only exists for perturbations 
of wavelength k = kq/A. The fastest growing mode, ei 
is expected to prevail in the case of a broadband initial 
fiuctuation spectrum. However it is possible to excite 
each mode independently, by imposing perturbations of 
appropriate wavelength. Excitation of ei triggers binary 
fusion of neighboring holes, leading to a final state with 
two holes shown in Fig. ll(a)-(c) which, owing to the 
periodicity of the system, is just the same final state as 
that obtained when N = 2. By contrast, excitation of 62 
leads to a more asymmetric hole-fusion senario, with sub- 
sequent nonlinear evolution leading to a one-hole state, 
see Fig. ll(d)-(e). 
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FIG. 8. (Color online) (a) TV = 2-cell system growth rate 
7 computed through Fourier-Hermite expansion with Nx = 
4:8, Va = 1.8, 8 = 0, 9 = —i and estimated from simulations 
with the Vlasov code VADOR. The red dashed lines represent 
the error bars in the growth rate estimated with VADOR. (b) 
Relative rate of change of 7 with A'',, . 



VI. CONCLUSIONS 

We have shown that the combination of spectral de- 
formation and Fourier-Hermite expansion can be used to 
compute the stability of nonlinear Vlasov-Poisson waves 
in an efficient and systematic manner. The computa- 
tion of unstable eigenfunctions for the BGK mode of 
Sec. V B illuminated the role of linear instability in trig- 
gering the subsequent vortex fusion. The ability to detect 
sub-dominant unstable modes permits a direct assesment 
of the role of perturbations of different wavelength in the 
evolution of the system. Our method will be used in a 
future paper to test the relevance of considering BGK 
equilibria in order to address Raman saturation. 

Spectral deformation was introduced to handle fila- 
mentation scale effects through damping of the advective 
term in Vlasov equation. This was the key to achiev- 
ing exponentially fast convergence of the eigenvalue com- 
putation and justify the choice of spectral methods. In 
contrast to estimates of the growth rate based on numer- 
ical integration with a Vlasov-Poisson solver, the com- 
bination of spectral deformation and Galerkin projection 
scales well with the dimension of phase space as the expo- 



nentially convergent series can be truncated very early, 
while still providing reliable results. We therefore find 
our method to be more suitable for extention to stability 
calculations in dimensions higher than 1-D, than an ap- 
proach based on direct integration with a Vlasov code. 
Moreover, there is in principle no obstacle to generaliz- 
ing the present method to the Vlasov-Maxwell system 
and to different geometries, as the choice of basis and 
transformed operator B{9) can be adjusted to the prob- 
lem at hand. Although our emphasis is on manipulating 
the structure of the linear operator A itself, rather than 
on a careful choice of basis, we expect that the method 
can be further optimized, if needed, through the tech- 
niques introduced in Ref. [21]. 



Appendix A: Hermite basis 

The Hermite polynomials H„ (u) we use here follow the 
standarization Hn{u) ^ 2"u" for large w, see Ref. [30]. 
Their explicit definition is 



H,{u) ^ (-l)^e" — e- . 

They satisfy the orthogonahty condition 

f +00 ^ 

H„,{u)Hn{u)e^''- du = (5,„„V7r2"n! , 

and therefore 



*'"(u)*„(u)du - S„ 



We note the following relations [30]: 

Hn+i{u) = 2uHniu) - 2nHn-iiu) . 



(Ala) 



(Alb) 



(A2) 



(A3a) 
(A3b) 



Using (A3) we can show that the orthonormalized 
Fourier-Hermite basis (12) satisfies 



(A4a) 



n + 1 n 
M*„(u) = W *„+i(u) + J -*„_i(u) . (A4b) 

Using the identity 

H,{u + c)=j2(i)Hk{u){2cy-', (A5) 
we can show 



fc=0 



+00 



du'^^u)^„{u + c) 



, j <n, 

1 , j ^n, 
§^Q)(-2c)^-", j>n. 

(A6) 
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FIG. 9. (Color online) Comparison of the isocontours of the 2-cell unstable eigenfunction as computed by (a) Fourier-Hermite 
method with N:, = 48, t)s = 1-8, iV„ = 109, 6 = (blue, dashed hue), iV„ = 110, 61 = (green, solid line), (b) Fourier-Hermite 
method with A^'^ = 48, Vs = 1.8, A^^ = 110, 9 = —\ (black, dashed line) and numerical integration of a small sinusoidal 
perturbation with the code VADOR with Ux = 1920 and Uv = 1000 (red, solid line), (c) The corresponding electric field 
eigenmode, using the same color-code as in Panels (a) and (b). 
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FIG. 10. (color online) (a) The unstable eigenmode of the 3- 
cell equilibrium, determined using = 72, = 110, «s = 
1.8, d = — i. (b) Predicted evolution in the linear neighbor- 
hood of the 3-cell equilibrium. 



Appendix B: Representation of translation operator 



The velocity translation operator U{9) is represented 
in Fourier-Hermite basis by the matrix elements 



A direct calculation gives, with the help of (A6), 

/ + 00 
-oo 

0, 
1, 



J <n, 
j = n, 



(B2) 



> n . 



However, in our computations we only need the ac- 
tion of U'^-'m„(6') on a vector /i*"" representing a function 
h{x, v). We can then exploit the sparsity of the represen- 
tation of the generator of velocity translations , 



id,, 



'j,m+l : 



and employ Krylov subspace approximations [26] to 
U{e)h{x,v) = e^^/i(x,w), to compute U'=^;„„(6l)/i™" in 
a stable and efRcient manner. 



Appendix C: Truncation 

In practice the infinite ladder of equations (15) has to 
be truncated by setting /^-+iJ(f) = /''■'""+i(t) = for 
some cutoff values and Ny . We present the truncated 
form of (21) to emphasize the treatment of boundary 
terms in Fourier space and to introduce the notations 
required to discuss the computation of the expansion co- 
efficients in Appendix D. For clarity, we only present the 
^ = case here; extension to the spectrally deformed 
case is straightforward. 

The truncated Fourier-Hermite expansions read 



f{Xj,Vm,t) 



A'x/2 



/"$,(xj)*,(u™), (Cla) 



dx / d«<i>'=(z)<i'J(«)c/(e„)$„(x)*„(«). 

-L/2 J-oo 

(Bl) 



E{xj,t)^ J2 E''<i>r{xj), (Clb) 
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FIG. 11. (color online) (a) The most unstable eigenmode ei of the 4-cell equilibrium, determined using A^^, — 96, Nv — 90, Vs = 
1.8, 6 = — i. (b) Linear evolution of a perturbation initialized along ei and (c) final state {t ~ 1600) reached after hole merging 
(using code VADOR). (d) The second unstable eigenmode 62 of the 4-cell equilibrium, (e) Linear evolution of a perturbation 
initialized along e2 and (f) final state {t ~ 1600) reached after hole merging (using code VADOR). 



j = 1, . . . , iV^ , m = 0, . . . , 7V„ , and = e''^'=«^, ko = 

2tt/L and a;„ — nL/Nx- The calculation of coefRcients 
is presented in Appendix D. 

Note that when taking odd derivatives of (Cla) with 
respect to Xj reahty of the result is not ensured, since the 
term — i(Af^/2) e~'(^="/^^'^''^j is not included in the sum. 
Thus, we have to set the (presumably small) Nx/2 term 
in the discrete Fourier transform of odd derivatives equal 
to zero, see Ref. [31, Chapter 3]. Then, proceeding as in 



Sec. IV with 6* = we get, when j > 1, 
while 



fcO 



Ini 



^ {-ikko8ki^5im, k^Nx/2, 
"jo, k = Nj2, 



and we have introduced 



(C2) 



(C3) 




0, k = N.J2 , 

-^n^J^ft'-" Wi> k ^ I, and |fc - /| < ^, 



otherwise , 



kn I 



Somfo , I ^OJ^ and 
otherwise . 



2 < k I < 2 



(C4) 

(C5) 
(C6) 



To store A^i im in a matrix Apj we substitute pairs 
of indices {k,j) and {l,m) with collective indices f3 — 
{k-l){N.^+l)+j and 7 = (/-l)(iV^, + l)+TO, respectively. 



Appendix D: Calculation of Fourier-Hermite 
coefficients 



Calculation of Fourier-Hermite coefficients f^'^ through 
a discrete analog of 

f^^ = i dx duf{x,VsU,t)<^>''{x)-i'^{u), (Dl) 

^0 J -00 
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where $'^(a;) = ig"''^'^''^, requires some care. Passage 
from an integral in a; to a discrete transform is readily 
handled by a standard Fast Fourier Transform (FFT). 
The analog of the discrete Fourier transform for inte- 
gration that involves an exponentially decaying function 
h{x) is Hermite-Gauss quadrature rule [30, 32] 

/ h{u)du ~ Wih{ui) , 

where the Gauss-Hermite weights are given by 



(iV + l)[vI/^(7.0] 

The sum has to be evaluated at an appropriate set of 
points Ui, the Gauss-Hermite quadrature points or ab- 
scissas. The quadrature points are given by the roots of 
HpiJ^i = and are not available in closed form. They 
can be evaluated as the solution to a symmetric, tridi- 
agonal eigenvalue problem and are not equispaced [33]. 
Evaluation of quadrature weights (D2) in a stable fashion 
is feasible through the use of Hermite recursion relations 
(A4b), rather than by direct use of (D2), see Refs. [32, 
Chapter 5] and [33, Remark 4.2]. 

The Fourier-Hermite coefRcients arc then given by 

AT, N 
n=l i=0 

where k = + 1, . . . , $'=(a;) = ]^e-"='=«"^, 

j = 0, . . . , N. In most cases the choice N = Ny provides 
a good approximation of Z*"'^, see Ref. [16, Chapter 4]. 
In some cases however, one might want to consider using 
N > Ny, for reasons related to some special features of 



the Hermite collocation grid. As already mentioned, the 
spacing between adjacent Hermite collocation points is 
not constant and thus the grid can by rarefied in regions 
where the distribution function has significant features, 
such as trapping areas. Moreover, the average grid spac- 
ing decreases only as 1/VN with increasing N, while at 
the same time the grid span increases as Vn (see [16, 
Chapter 17]). Therefore, we have to increase N signifi- 
cantly to ensure adequate resolution in areas of interest, 
while at the same time we are forced to include contri- 
butions from many large v collocation points, where the 
distribution function is practically zero. 

We overcome this problem by a non-standard scheme 
for the computation of f'^^. Distribution functions of 
interest fall off faster than any polynomial in |u| for 
large v. We can therefore introduce appropriate veloc- 
ity cutoff values u™'" < < u"^^^ and replace X)iLo 
(Dl), with E"=«i' 'Where m"''^ < u„, and < u™^'^. 
One can determine ni, n-i by counting eigenvalues within 
of the eigenvalue problem to which Hat+i = 
reduces, as prescribed in Ref. [26, Lecture 30]. We can 
now pick N 3> iV^, ensuring adequate resolution within 
the region of interest, without unnecessarily increasing 
the truncation order Ny. One still has to solve the sym- 
metric, tridiagonal eigenvalue problem for the colloca- 
tion points of order N . However, we only need a small 
subset u™'" < u < u™^'' of the spectrum and an effi- 
cient algorithm based on bisection can be employed [34] . 
Appropriate m™™, u^'^^ and TV have to be determined 
empirically for the problem at hand. Good starting val- 
ues for u™'" and n^'^^ are uq and ujv„ respectively, since 
Hermite polynomials of order Ny do not oscillate outside 
(uo,WArJ. 

To transform from Hermite coefficient space back to 
physical space efficiently, Clenshaw recurrence [32, Chap- 
ter 5] can be used following the FFT in (Cla). 
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